Objetivo de aprendizaje
Este caso práctico muestra como leer y graficar datos espaciales en R. Para ello, vamos a trabajar con los datos de temperatura mínima registradas en España por las estaciones metereológicas de la AEMET.
Tarea 1: Abrimos RStudio
El presente análisis se va a realizar empleando RStudio, por lo que empezaremos abriendo el programa y creando un nuevo script de R en File>New File>R script.
Tarea 2: Importamos y describimos los datos objeto de estudio
El conjunto de datos proporcionado (tempmin.csv) contiene el nivel de
temperatura del aire en España entre el 6 y el 10 de Enero de 20211. Estos
datos han sido descargados usando la librería climaemet (Pizarro, Hernangómez, & Fernández-Avilés, 2021) y han
sido posteriormente tratados para su uso en esta práctica.
El primer paso consiste en importar la base de datos de temperatura mínima. El
archivo está en formato csv, por lo cual es un fichero de texto plano. Podemos
usar varias funciones para realizar la importación, en este caso vamos a emplear
paquetes del tidyverse para realizar todo el tratamiento de datos:
library(readr)
tmin <- read_csv("data/tempmin.csv")
Q1: ¿Qué información tiene tmin?
| fecha | indicativo | tmin | longitud | latitud |
|---|---|---|---|---|
| 2021-01-06 | 4358X | -4.7 | -5.880556 | 38.95556 |
| 2021-01-06 | 4220X | -7.0 | -4.616389 | 39.08861 |
| 2021-01-06 | 6106X | 4.7 | -4.748333 | 37.02944 |
| 2021-01-06 | 9698U | -6.8 | 0.865278 | 42.20528 |
| 2021-01-06 | 4410X | -3.4 | -6.385556 | 38.91583 |
| 2021-01-06 | 1331A | 1.0 | -7.031389 | 43.52472 |
Podemos observar que la tabla generada contiene 5 columnas distintas:
fecha: Indicando la fecha de observación.
indicativo: Es el identificador de la estación de la AEMET que registró el
dato.
tmin: Dato de temperatura mínima registrada en cada fecha por la estación
correspondiente.
longitud,latitud: Coordenadas geográficas de la estación
Tarea 3. Convertir tmin a un objeto de la clase espacial geoR
Para convertir un objeto a geodata (el formato requerido por geoR),
proporcionaremos una tabla con las coordenadas y los valores a incluir en cada
coordenada. En este ejemplo, vamos a emplear sólamente los datos
correspondientes al 8 de enero.
library(dplyr)
library(geoR)
tmin_geoR <- tmin %>%
filter(fecha == "2021-01-08") %>%
# Seleccionamos las columnas de interés
dplyr::select(longitud, latitud, tmin) %>%
# Y creamos el objeto geodata
as.geodata(
coords.col = 1:2,
data.col = 3
)
summary(tmin_geoR)
#> Number of data points: 211
#>
#> Coordinates summary
#> longitud latitud
#> min -9.291389 35.27639
#> max 4.215556 43.78611
#>
#> Distance summary
#> min max
#> 0.01024389 13.85144264
#>
#> Data summary
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -14.9000000 -4.6000000 -0.5000000 -0.6293839 3.5000000 13.6000000
plot(tmin_geoR)
Figure 1: Objetos en geoR
DIEGO, y yo por qué se que España es esto 4326…
Está explicado en la sección de CRS, es el código del CRS que usan los GPS, el habitual para coordenadas longitud/latitud
Tarea 4. Convertir tmin a un objeto de la clase espacial sf
En esta tarea, convertiremos los datos de tmin en un objeto espacial sf, es
decir, datos espaciales de tipo vector.
Los datos de tmin contienen coordenadas geográficas longitud/latitud, asi que
como vimos en la sección Sistema de Referencia de Coordenadas (CRS) el CRS a
emplear ha de ser un CRS geográfico. Usaremos el código EPSG 4326, que
corresponde a coordenadas geográficas y suele ser el habitual en este tipo de
situaciones.
library(sf)
tmin_sf <- st_as_sf(tmin,
coords = c("longitud", "latitud"),
crs = 4326
)
tmin_sf
#> Simple feature collection with 1066 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -9.291389 ymin: 35.27639 xmax: 4.215556 ymax: 43.78611
#> Geodetic CRS: WGS 84
#> # A tibble: 1,066 x 4
#> fecha indicativo tmin geometry
#> * <date> <chr> <dbl> <POINT [°]>
#> 1 2021-01-06 4358X -4.7 (-5.880556 38.95556)
#> 2 2021-01-06 4220X -7 (-4.616389 39.08861)
#> 3 2021-01-06 6106X 4.7 (-4.748333 37.02944)
#> 4 2021-01-06 9698U -6.8 (0.865278 42.20528)
#> 5 2021-01-06 4410X -3.4 (-6.385556 38.91583)
#> 6 2021-01-06 1331A 1 (-7.031389 43.52472)
#> 7 2021-01-06 1690A -0.1 (-7.859722 42.32528)
#> 8 2021-01-06 8489X -8 (-0.255833 40.43333)
#> 9 2021-01-06 8025 2 (-0.494167 38.3725)
#> 10 2021-01-06 9784P -10 (0.224722 42.63)
#> # ... with 1,056 more rows
Tarea 5: Dibujemos las estaciones de monitoreo de la temperaria mínima en un mapa de España. Ámbito espacial.
Vamos además a incluir una capa de las comunidades autónomas de España. Para
ello usaremos un paquete API que nos proporciona esta información en formato
sf:
library(mapSpain)
# sf object
esp <- esp_get_ccaa() %>%
# No vamos a usar Canarias en este análisis
filter(ine.ccaa.name != "Canarias")
plot(esp$geometry) # Dibujamo el mapa de España menos las Islas Canarias
Figure 2: Mapa de España (Sin Canarias)
Q2: ¿Tengo el Sistema de referencia de coordenadas (CRS) de las estaciones de monitoreo en la misma proyección que el contorno de España?
Como se comentó en la sección correspondiente, cuando se emplean datos geográficos provenientes de varias fuentes, es necesario asegurarse de que ambos objetos están usando el mismo CRS. Vamos a comprobarlo:
st_crs(tmin_sf) == st_crs(esp)
#> [1] FALSE
Vemos que no lo están, por lo que vamos a proyectar las coordenadas a un CRS
común. En este caso usaremos el CRS de referencia de tmin_sf:
esp2 <- st_transform(esp, st_crs(tmin_sf))
st_crs(tmin_sf) == st_crs(esp2)
#> [1] TRUE
Dibujamos las estaciones de monitoreo con el contorno de España. Vamos a usar el
paquete ggplot2 como referencia, sin embargo existen varios paquetes
especializados en mapas temáticos, como pueden ser tmapo mapsf.
library(ggplot2)
ggplot(esp2) +
# Para graficar objetos sf debemos usar geom_sf()
geom_sf() +
geom_sf(data = tmin_sf) +
theme_light() +
labs(
title = "Estaciones de monitoreo AEMET en España",
subtitle = "excluyendo las Islas Canarias"
) +
theme(
plot.title = element_text(
size = 12,
face = "bold"
),
plot.subtitle = element_text(
size = 8,
face = "italic"
)
)
Figure 3: Estaciones de AEMET en la Península Ibérica
Q3. Mis datos y el contorno de España están en el mismo CRS, pero ¿es el adecuado?
DIEGO, POR QUÉ TRANFSORMAMOS, PARA PASAR A UTM E INTERPRETAR EN METROS?
Explicado más adelante, para hacer el grid en metros cuadrados. Si no, se haría en grados cuadrados (¿?) porque long/lat no es medida de superficie
Tarea 6: Representamos la variable temperatura mínima tmin para el día 8 de
enero de 2021.
En la siguiente tarea, seleccionaremos los datos correspondientes al 8 de enero de 2021 y crearemos un mapa temático en el que representaremos los valores de temperatura mínima registrados en cada estación mediante un código de colores.
tmin_8enero <- tmin_sf %>%
# seleccionamos el día y la variable
filter(fecha == "2021-01-08")
plot(tmin_8enero["tmin"],
main = "Temperatura mínima (8-enero-2021)",
pch = 8
)
Figure 4: Mapa de puntos con temperatura mínima
Podemos utilizar el ámbito espacial, el contorno de España para graficar y contar la historia de la Filomena un poco mejor.
Figure 5: Mapa completo con temperatura mínima
Q4. El mapa ha quedado muy claro. Vemos como los datos nos cuentan la historia de Filomena en aquellos sitios donde se tomaron mediciones, pero ¿podríamos tener un mapa de interpolación para tener una estimación de la temperatura mínima en las partes donde la AEMET no tiene estación de monitoreo?
Tal y como se avanzó en teoría, parece lógico pensar que aquellos puntos que estén cerca tendrán valores similares así que tomemos ventaja de la dependencia espacial y utilicemos un método determinista, como la Distancia Inversa Ponderada, comúnmente conocido por su acrónimo inglés IDW (Inverse distance weighted), el cual es uno de los métodos más simples para llevar para llevar a cabo una interpolación espacial.
En este tipo de análisis, es crucial que el CRS sea el apropiado. En este caso, ya definimos el CRS como un CRS geográfico (es decir, usando coordenadas de longitud y latitud). Sin embargo, para el ejercicio de interpolación es más adecuado usar un CRS local (que provoque pocas deformaciones en la proyección de España) y en alguna unidad de distancia, como metros (ya vimos que en los CRS geográficos las unidades son grados).
Si usamos el paquete crsuggest podemos observar los CRS sugeridos:
library(crsuggest)
sugiere <- suggest_crs(tmin_8enero, units = "m", limit = 5)
# Usamos la sugerencia del paquete
crs_sugerido <- st_crs(sugiere[1, ]$crs_proj4)
esp3 <- st_transform(esp2, crs_sugerido)
tmin_8enero3 <- st_transform(tmin_8enero, crs_sugerido)
A continuación, para realizar la interpolación, necesitamos generar una malla que representará las celdas de las que queremos obtener el valor interpolado.
Dado que hemos proyectado nuestros datos a un CRS cuya unidad son los metros, podemos definir el tamaño de cada celda en metros cuadrados. En este caso vamos a usar celdas de 49 kms cuadrados (7 x 7 kms):
malla_sf <- st_make_grid(
esp3,
cellsize = 7000
)
Graficamos la superficie para ver exactamente qué hemos construido:
ggplot(esp3) +
geom_sf() +
geom_sf(
data = malla_sf,
size = 0.1,
col = "red", alpha = 1,
fill = NA
) +
geom_sf(
data = tmin_8enero3,
aes(fill = "AEMET Stations"), size = 5, shape = 21,
color = "blue"
) +
scale_fill_manual(values = adjustcolor("blue", alpha.f = 0.2)) +
theme_void() +
theme(legend.position = "bottom") +
labs(
title = "Cuadrícula espacial para interpolar",
fill = ""
)
Figure 6: Malla de puntos para interpolación
Se puede observar claramente cada una de las celdas que hemos creado. La interpolación asignará un valor a cada uno de ellas.
A continuación podemos llevar a cabo la interpolación usando el paquete gstat.
Además, en lugar de celdas (polígonos) es necesario usar puntos en la
interpolación. Calcularemos por tanto un punto representativo de cada celda, el
centroide, que es el punto resultante de realizar la media arimética de las
coordenadas de loss puntos que componen los lados de cada celda
# Calculamos centroide
malla_sf_cent <- st_centroid(malla_sf, of_largest_polygon = TRUE)
library(gstat)
tmin_idw <- idw(
# Indicamos la variable que queremos interpolar
tmin ~ 1,
# Indicamos el conjunto de datos donde está la variable
tmin_8enero3,
# Indicamos la malla de destino, en sf
newdata = malla_sf_cent,
idp = 2.0 # Especifica la potencia de la IDW
)
#> [inverse distance weighted interpolation]
head(tmin_idw)
#> Simple feature collection with 6 features and 2 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 145790.9 ymin: 68957.31 xmax: 180790.9 ymax: 68957.31
#> CRS: +proj=lcc +lat_1=40 +lat_0=40 +lon_0=0 +k_0=0.9988085293 +x_0=600000 +y_0=600000 +a=6378298.3 +rf=294.73 +pm=madrid +units=m +no_defs
#> var1.pred var1.var geometry
#> 1 2.513414 NA POINT (145790.9 68957.31)
#> 2 2.572912 NA POINT (152790.9 68957.31)
#> 3 2.633632 NA POINT (159790.9 68957.31)
#> 4 2.695589 NA POINT (166790.9 68957.31)
#> 5 2.758802 NA POINT (173790.9 68957.31)
#> 6 2.823285 NA POINT (180790.9 68957.31)
Representamos la interpolación con un mapa y mapa de contorno muy utilizado para
representar datos espaciales. Para ello, vamos a usar el paquete raster
convirtiendo nuestro objeto interpolado.
# Convertimos de sf a SpatiaPixels
# Esto funciona porque nuestros puntos sf están espaciados regularmente
tmin_pixels <- tmin_idw %>%
as("Spatial") %>%
as("SpatialPixels")
library(raster)
# Creamos un raster de nuestros pixels
rast_esp <- raster(tmin_pixels)
# Transferimos valores del objeto sf al raster
rast_esp2 <- rasterize(
tmin_idw,
rast_esp,
field = "var1.pred",
fun = mean
)
# Además, podemos recortar el raster a la forma de España
rast_esp_mask <- mask(rast_esp2, esp3)
plot(rast_esp_mask, col = colores)
contour(rast_esp2, add = TRUE)
Figure 7: Mapa raster con lineas de nivel
Podemos realizar el mismo mapa usando ggplot2 y la función
geom_contour_filled:
# Creo una tabla para geom contour
coordenadas <- st_coordinates(tmin_idw)
valor <- tmin_idw$var1.pred
idw_df <- data.frame(
# Necesitamos redondear las coordenadas
latitud = round(coordenadas[, 2], 6),
longitud = round(coordenadas[, 1], 6),
tmin = valor
)
ggplot() +
geom_contour_filled(
data = idw_df,
aes(x = longitud, y = latitud, z = tmin),
na.rm = TRUE,
breaks = cortes
) +
# Reajustamos la escala de colores
scale_fill_manual(values = colores) +
# CCAA
geom_sf(data = esp3, fill = NA) +
theme_minimal() +
theme(axis.title = element_blank()) +
labs(
fill = "Temp. (º)",
title = "Temperatura mínima interpolada",
subtitle = "8 de Enero 2021",
caption = "Datos: AEMET"
)
Figure 8: Mapa con ggplot2 y lineas de nivel
Representamos en 3D el mapa anterior, muy utilizado también en datos espaciales
DIEGO, se útiliza mucho en espacial pero no se si este caso es muy ilustrativo….r
Por mi lo quitamos
persp(rast_esp_mask, theta = 2000, d = 0.1)
Las fechas seleccionadas coinciden con el periodo en el que la tormenta Filomena tuvo su auge en la Península Ibérica.↩︎